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ABSTRACT 

cn : 

We investigate how galactic disks react to external tidal torques. We cal- 
culate the strength and radial dependence of torques on disks that arise from 
a misalignment between the disk and the main axis system of a flattened dark 

^ matter halo. Density profile, misalignment and flattening of the halo are chosen 

to match the corresponding values typically found for dark matter halos in large 
cosmological N-body simulations. We find that except for in the very inner re- 
\ gions, the torques are well- described by a power law of the form r oc r~ 2 - 5 . For 

torques as they arise in typical cosmological settings, the magnitude of the torque 
is large enough for the entire disk to react to the torque in less than the Hubble 
time. We demonstrate analytically that disks which are originally located in the 
xy-plane and which are subjected to a torque around the x-axis tilt around the 
y-axis, as also found in fully non-linear N-body simulations. We further demon- 
strate that that the torque causes the radius of a chosen particle to increase with 
time. Investigations of tilting disks which treat the disk as a set of solid rings 
thus may systematically overestimate the effects of the torque by a factor of two. 
For torques of the form we investigate, the inner regions of the disk react to the 
torque faster than the outer regions, resulting in a trailing warp. We then study 

. \ the effect of the self-gravity of the disk in such a scenario using numerical N-body 

models. Self-gravity flattens out the inner regions of the disk, but these regions 
are tilted with respect to their initial plane followed by a non-flat outer region 
whose tilt decreases with radius. The "warp radius," which marks the end of 
the inner flat disk, grows throughout the disk at a rate that depends only on the 
strength of the torque and the local surface density of the disk. 
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Subject headings: stellar dynamics — methods: N-body simulations — galaxies: 
evolution — galaxies: kinematics and dynamics — galaxies: spiral — galaxies: 
halos 

1. Introduction 

Many edge-on disk galaxies show integral-sign or S-shaped warps, where the majority 
of the disk is planar but where the outer region of the disk lies above that plane on one side 
of the galaxy and below the plane on the other (Binney 1992). The Milky Way is warped 
both in neutral hydrogen (e.g. Diplas & Savage 1991) and in the stellar distribution (Reed 
1996; Lopez-Corredoira et al 2002b). Reshetnikov & Combes (1998) estimate that half of all 
disk galaxies have optical warps, and most HI disks which extend beyond the optical disks 
appear to be warped (Bosma 1981; Briggs 1990; Christodoulou et al 1993). 

A number of methods have been proposed for creating and maintaining warps. Several 
authors (Lynden-Bell 1965; Sparke & Casertano 1988) have suggested that the system of 
discrete particles which make up the disk may have normal bending modes that could be 
excited. Battaner et al (1990) have investigated magnetic fields as a cause of warps, while 
Toomre (1983) and Dekel & Shlosman (1983) suggested that disks askew in flattened dark 
matter halos could develop long-lived warps, assuming the radial profile of the halo were 
appropriately fine-tuned. 

Following on the idea that infalling material will shift the angular momentum of a galaxy 
(Ryden 1988; Quinn & Binney 1992), Ostriker & Binney (1989) studied how a disk of massive 
rings reacts to a slewing disk potential. They found that the reaction of the disk depends 
largely on its surface density, with warps appearing in regions of low surface density. Also 
motivated by the cosmic infall of material with differing angular momentum, Debattista & 
Sellwood (1999) found that when the angular momentum of a halo and disk are misaligned, 
dynamical friction transfers angular momentum between them in such a way as to produce 
a long-lived warp. 

More recently, Lopez-Corredoira et al (2002a) examined the torques produced by the 
transfer of angular momentum from gas falling onto a galactic disk, and the warping of the 
disk in response to these torques as well as the internal torques due to the interaction of 
different rings within the disk. They find a disk response which mirrors the Milky Way warp, 
though the effects of embedding the disk in a dark halo are as yet insufficiently modelled. 
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In a cosmological setting, a galactic disk experiences gravitational torques from three 
different sources: external galaxies, non-spherically-symmetric substructure in the halo, and 
misalignment between the disk and halo; if the dark matter halo is not spherical, as suggested 
both by simulations (e.g. Frenk et al 1985; Katz 1991; Dubinski & Carlberg 1991; Cole & 
Lacey 1996) and observations (Sackett 1999), then it will exert a gravitational torque on a 
disk not in its symmetry plane. 

In this paper, we evaluate the effect of the gravitational tidal torques a typical galactic 
disk experiences from its misalignment with the halo, and study whether these torques 
provide a possible origin for warped disks. We first develop a framework for how orbits 
in a disk react to torques. We then investigate the form and magnitude of torques due to 
flattened misaligned halos. The effect of these torques on massless disks is examined by 
analytic and numerical techniques. Finally, we perform N-body simulations of massive disks 
resembling the Milky Way and investigate how the self-gravity of the disk may affect its 
reaction to cosmological torques. 



2. Origin of Torques on Galactic Disks 
2.1. Framework for Torqued Orbits 

2.1.1. Introduction to torque 

The detailed reaction of a stellar disk to a torque is derived in section 3. Here we 
introduce a simplified model to establish a framework in which the reaction of the disks to 
the torques can be studied, and define the terms and symbols we use. Consider a star in the 
disk on a circular orbit in the x|/-plane around the origin with angular velocity u(r) = v c {r)/r 
(see Figure 1). We apply a torque around the a;-axis by accelerating the star in the y and z 
directions around the a;- axis with an antisymmetric pseudo-tensor Ty. If the acceleration in 
the i direction is and the j-position of the star is r,j then the applied acceleration is 

ai = J2 T ii r j (!) 

j 

where 

/ \ 
T = -r(r) . 
\ r(r) / 

This is an angular acceleration (or equivalently a specific torque) around the x-axis of magni- 
tude a/r = r(r). Since the angular acceleration is the time derivative of the angular velocity, 
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the x-component of the angular velocity increases in response to the torque, and the angular 
momentum of the star will slew toward the x axis, tilting the orbit toward the x-axis. 



2.1.2. Direction of tilt 

It is worth taking a minute here to elaborate on the point that the orbit tilts toward 
the x axis. One might think that a torque about the x axis would tilt the plane of the orbit 
around the x axis, rather than toward it. Indeed, this is what happens initially, and seems 
to be the thinking behind attempts to explain the Milky Way warp using the Magellanic 
Clouds (e.g. Weinberg 1998; Tsuchiya 2002). However, Garcia- Ruiz et al (2000) examined 
the Lagrangian of the system and conducted N-body experiments to demonstrate that a 
warp caused by the Magellanic Clouds would be perpendicular to the observed Milky Way 
warp. 

This result can be understood by examining the equations of motion of the system for 
a particle on a circular orbit in the xy plane about the z-axis: 

^!-!^ r + T xr (2) 

dt 2 r dr 

If the untorqued orbit has an angular frequency u, and the torque is directed around the 
x-axis, then the equations of motion for the three Cartesian coordinates are 

-uo 2 x (3) 



dt 2 
fy 
dt 2 

cPz 

dt 2 



—u) 2 y — tz (4) 
-uj 2 z + ry. (5) 



The x motion is decoupled from the other coordinates, and therefore is unaffected by the 
torque. By substituting equation (4) into equation (5), the system of equations can be 
converted into a single fourth-order ordinary differential equation for z of the form 

^ + 2u ^ + (^ + ^ = 0. (6) 

Using the ansatz z = Ae im , the roots of the characteristic equation are 

n 2 = u 2 ± ir. (7) 

The two solutions correspond to growing (— ) and decaying (+) modes. The growing mode 
will dominate over time, so we neglect the decaying mode. As d 2 z/dt 2 = — f2 2 z, equations (5) 
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Fig. 1. — Two views of a circular orbit with angular velocity vector on tilted from the xy- 
plane toward the a;- axis by an angle 9. is the azimuthal angle around the orbit. A torque 
r directed along the x-axis produces forces F, as in equation (1). 



Acceleration Orbit 




y 



Fig. 2. — Left: Acceleration vectors due to the torque for an orbit in the xy-pl&ne. Right: 
The z response lags y by one quarter of an orbit, causing the orbit to tilt around the y axis 
toward the x axis. 
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and (7) yield 

V = iz. (8) 

For an oscillating y, this means that the phase of y is a quarter of an orbit ahead of z, or 
that the phase of z is a quarter of an orbit behind y. As demonstrated in Figure 2, z reaches 
its maximum a quarter of an orbit after y, causing the orbit to tilt around the y-axis toward 
the a;-axis. 

Therefore, we expect that any tilting (and therefore warping, which is a tilt that varies 
with radius) will occur toward the plane perpendicular to the torque. Since a torque is a 
transfer of angular momentum, a simple way to understand the result of Garcia-Ruiz et al 
(2000) is that the galactic disk acquires angular momentum with the same direction as the 
orbital angular momentum of the satellite, and therefore it will warp toward the plane of 
the satellite's orbit, not toward the satellite itself. 



2.1.3. Tilting timescale 

For an orbit in the xy-plane, the angular momentum is initially aligned with the z axis 
(w = u (r) z). The torque adds angular momentum in the x direction at a constant rate: 

du> 

— — = T(r) x. 
dt v 1 

If the radius of the orbit stays constant, then the angular velocity grows as 

u(t) = tr(r) x + Ld (r) z. 

As the torque adds angular momentum along the x axis, the the orbit tilts toward the 
x axis by an angle 9, where 

uj x tr(r) ,„ x 
tan# = — = —7-4- 9 

lo z LO Q {r) 

The rate of tilting at small angles is 

d9 Hr) 



dt oj${r) 

with a characteristic timescale to tilt one radian of 



<„■« = do) 

r(r) 



Another motivation for this form is that uj is the angular momentum per unit mass 
and r is the specific torque, which is the rate of change of the specific angular momentum. 
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Note that even for a flat rotation curve where uo = v c /r, t t iit is a monotonically increasing 
function of r if r(r) falls off faster than r -1 ; i.e. for a torque that decreases with radius, 
the inner regions of the disk should always tilt faster than the outer regions, resulting in a 
trailing warp. 



2.2. Expected torques from misaligned flattened halos 

There is observational (Sackett 1999) and theoretical (Katz 1991; Cole & Lacey 1996) 
evidence that dark matter halos are flattened with c/a axis ratios ranging as low as 0.5. 
Cosmological simulations (Frenk et al 1985; Katz 1991; Warren et al 1992; Porciani et al 
2001) show that while the angular momentum and minor axis of a flattened halo will be 
correlated, there will be significant misalignments between them in a large number of halos. 
Recent high-resolution simulations suggest that the angular momentum of the baryons may 
not even be aligned with that of the dark matter (van den Bosch et al 2002), and that 
therefore the disks that form from those baryons are usually misaligned with the minor axis 
of the halo matter distribution. 

We calculate the magnitude of the torque for a disk misaligned in a flattened halo with 
a radial dependence of an NFW profile (Navarro et al 1996) but flattened along the z axis: 

p(x ' y ' z) = ( i w7T / ^ ( n ) 

\m/r s ) (1 + m/r s ) 

for a modified radius m 2 = x 2 + y 2 + z 2 /q 2 where q is the c/a axis ratio. The force from 
this distribution on a given point is calculated using equation (2-88) of Binney & Tremaine 
(1987). For reference, we also calculate the torques inside flattened isothermal profiles, as in 
equation (2-54a) of Binney & Tremaine (1987). The torques are calculated by evaluating the 
forces at opposite points along a fictitious disk centered at the origin and placed at an angle 
9 to the xy plane (see Figure 1). The acceleration at radius r due to a torque of magnitude 
r is out of the plane, is of opposite sign on opposite sides of the disk, and is of magnitude 
F = rr. Given the forces on two test points F\ and F 2 at opposite sides of the disk, we take 
the component of the force out of the plane of the disk F ±1 and F ±2 . The torque is 

F±i — F ±2 

The angular velocity to is determined at each point from the radial force: uo = a/ F r /r where 
the radial force F r is defined to be positive when it is directed toward the center of the halo. 
The average value of uo for the two points across the disk is used. The tilting timescale is 
given by equation (10). 
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Fig. 3. — The symbols indicate the torques experienced by disks misaligned inside flattened 
halo profiles, expressed in terms of the tilting timescale t t iit- The lines are power law fits 
for log/? > 0.1. Note that stronger torques have shorter timescales and appear lower in 
these graphs. Top-left: Torques inside NFW halos with virial velocities t> 20 o = 175 km s -1 , 
flattenings q = 0.7, and concentration parameters c 2 oo — 15, for disks misaligned by 10°, 20°, 
and 30° from the symmetry plane. Top-right: Torques inside NFW halos with virial velocities 
^200 — 100, 175, and 200 km s -1 , flattenings q = 0.7, and concentration parameters c 2 oo = 15, 
for disks misaligned by 20° from the symmetry plane. Bottom-left: Torques inside NFW halos 
with virial velocities t> 2 oo = 175 km s -1 , flattenings q = 0.5, 0.7, and 0.9, and concentration 
parameters c 20 o — 15, for disks misaligned by 20° from the symmetry plane. Bottom-right: 
Torques inside isothermal halos with virial velocities t> 20 o = 175 km s -1 , flattenings q = 
0.5, 0.7, and 0.9, and core radii R c — 1 kpc, for disks misaligned by 20° from the symmetry 
plane. 
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Figure 3 illustrates the torques derived using this method for halos with a cross-section 
of properties expected for galactic halos. The fiducial NFW model has flattening q = 0.7, 
virial velocity v c = 175kms _1 (corresponding to v rot = 205 km s" 1 at a galactocentric radius 
of 10 kpc), and a disk at an angle 20° to the symmetry plane of the halo. All NFW halos 
have concentration parameters C200 — 15 (Navarro et al 1997). The different panels show the 
effect of varying the disk angle by 10°, the axis ratio from 0.5 to 0.9, and the virial velocity 
by 75 km s^ 1 compared to the fiducial model. The bottom-right panel of Figure 3 shows 
the torques inside an isothermal profile with a range of axial ratios. The magnitude of the 
torque in the isothermal case is similar to that in the NFW case, but falls off less rapidly 
with radius. 

Figure 3 demonstrates that the torque from a misaligned flattened halo is significant 
to its evolution on a cosmological timescale, as the tilting timescale i t iit is l ess then the 
Hubble time over the entire length of the disk for all halos modelled. The torques scale 
approximately as the density, and therefore fall off with radius as r(r) oc r 13 with (3 ranging 
from —1 at small radii to —3 at large radii for NFW halos and equal to —2 for isothermal 
halos. The typical torque follows the relation 



The angular velocity iv(r) oc r a where a ranges from —0.5 to —2 in NFW halos, and is —1 in 
isothermal halos. The tilting timescale t t iit, as calculated using equation (10), rises almost 
linearly with radius; the timescale is shorter in the inner regions of the disk, and therefore 
the inner disk tilts faster than the outer disk in response to torques from misaligned halos. 
For example, in an isothermal halo r(r) oc r~ 2 while u(r) oc r _1 so t t nt = u(r)/r(r) oc r. 
The torque increases more rapidly toward the center than does the disk's ability to resist 
the torque due to its angular momentum. 

The further the halos and disks are misaligned, the stronger the torque is; however, 
the torque profile is mostly unchanged for angles beyond 20°. Increasing the virial velocity, 
and therefore the mass of the halo, increases the magnitude of the torque. Because the 
concentration C200 is held constant amongst these models, changing the virial velocity also 
changes the scale radius r s , which can be seen in the shifting of the radius of the knee in 
the profiles of the top-right panel of Figure 3. The flattening of the halo has a large effect 
on the magnitude of the torque, with torques strengthening as the halo departs further from 
spherical symmetry. The torques fall off more slowly in the isothermal halos than in the 
NFW halos, but are of similar magnitude over most of the disk radius. 




(12) 
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3. Orbits in Torqued Disks 
3.1. Reaction of a massless disk 



The results of section 2.1 assume that the torque does not change the radius of the orbit. 
This is a poor approximation, since the torque does work on the star whenever the torque 
is not orthogonal to the orbital angular momentum. In this section we do a more thorough 
analysis, taking into account the changing radius of the orbit. 

We define 9 as in section 2.1, and as the azimuthal angle around the orbit (see 
Figure 1). Based on the results of section 2.2, the torque can be expressed as a power law: 

r = (13) 

which exerts a force 

F = mr x r 

( \/3+i (14) 
= mr r ( ^ J (cos 4> sin 9 y + sin (f) z) . 

This force torques the orbit and alters the angular momentum L. If the orbit can be 
considered to be a circular orbit with circular velocity v c (r) (where the radial dependence is 
a function of the underlying potential), then 



L = r x p = mrv c (r) u) (15) 

where cD is the unit vector in the direction perpendicular to the orbit. L has an additional 
component of order mr 2 {d9 j dt)y because the orbit is tilting, but it does not provide an 
independent constraint, so we neglect it. The time derivative of L is 



§ = [mv c (r) sin 9% + mr sin 6%% + mrv c (r)f t cos 9} x 
+ [mv c (r) cos 9^ + mr cos 9^^ - mrv c (r)^ sin 9] z. 

The change in angular momentum comes from the torque N imparted to the system: 
N = r X F. Taking F from equation (14) and averaging over one orbit to remove the 
dependence (assuming that dt oc d(f> over one orbit, i.e. that cu(r) is slowly varying so the 
0-averaged torque is equal to the time-averaged torque) gives 

(AO = \mr rl U) (1 + sin 2 9) x 

( \/3+2 U') 

+ \mT§r\ sin 9 cos 9 f ^ J z 
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Setting dL/dt = (N) and solving for dr/dt and d9/dt, 

dr 



Tt = T ° r »UJ 



r \ 2 \ dvc 1 ! 

v c (r) + r— 
dr 



sin 9 



d6 
~dt 



1 r r 



2v c (r) \r J 



— COS0. 



(18) 
(19) 



We could not find an analytic solution to this coupled pair of differential equations, so 
they need to be solved numerically. However, if dr/dt is small then equation (19) can be 
solved at a fixed radius. The solution is 



2 cos 



6(r,t) 

where the timescale at a given radius is 

to(r) = 



1 1 + e 



t/t (r) 



V2 Vl + e 2 */*oW 



(20) 



2v c (r) 
Wo 



-(3-1 



= 2 M r ) 

r(r) 

and the circular velocity v c (r) depends on the halo potential. 



(21) 



(22) 



Comparison of equations (22) and (10) reveals that the main effect of allowing a second 
degree of freedom, by permitting r to vary, is to increase the tilting timescale by a factor 
of two. This factor is also evident from a comparison of the dashed and the solid line in 
Figure 4. Theoretical studies of warps which artificially restrict this degree of freedom, for 
example by approximating the disk of stars by solid rings of constant radius, systematically 
overestimate the effect of the torques on the disk by a factor of two. 



3.2. Comparison with simulations 

The simplified model of equation (9), the detailed derivation (20), and the numerical 
solution to equations (18) and (19), provide three predictions of how a torqued orbit should 
evolve. To see how well these predictions work, we compare them to a simulated disk of 
massless particles subject to a torque that takes all the forces explicitly into account. 

We perform a simulation of a disk of 16000 massless tracer particles placed in 800 con- 
centric rings at radii spaced 0.05 kpc apart, each containing 20 particles. They are evolved for 
1 Gyr using a predictor-corrector code under the influence of an NFW potential (c 2 oo = 15, 
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R [kpc] 

Fig. 4. — Tilt versus radius for rings in a disk subject to a torque of the form (13) with 
(3 = —2.5 and r = 10~ 30 s -2 at r = 1 kpc for 1 Gyr. Filled diamonds correspond to a 
N-body simulation (binned in spherical shells 3 kpc wide), open squares to the numerical 
solution of equations (18) and (19), the dashed line to the analytic prediction from the 
simplified model of equation (9), and the solid line to the analytic prediction of equation (20). 
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r 2 oo = 200 kpc) and a torquing force of the form (12). We cap the torque at r = t in- 
side r in order to prevent extremely small timesteps for particles near the center where the 
torquing force diverges. The particles were then binned into spherical shells 3 kpc wide. We 
calculated the moment of inertia tensor 1^ = Xlfc m fc x «,fc x j,fc for the particles within each 
shell and diagonalized it to find the direction of the minor axis (the axis with the smallest 
moment of inertia). The minor axis was found to be tilted toward the x-axis by an angle 9 
that varied from bin to bin. 

We also evolved the coupled set of differential equations (18) and (19) forward 1 Gyr 
with the same values of f3 and To(ro) as above using a fourth-order Runge-Kutta integrator 
for a set of initial conditions 9 = 0, r = 0, 0.5, 1, 1.5 . . . 100 kpc representing an initially flat 
disk. 

Figure 4 compares the simulation (filled diamonds) with the numerical solution to equa- 
tions (18) and (19) (open squares), the simplified model (9) (dashed line), and the detailed 
analytic prediction (20) (solid line). It is apparent that at most radii, (20) is a good ap- 
proximation to the simulation — much better than (9), demonstrating the importance of 
including the radial drift of the orbits. At small tilt angles, the difference between the solid 
line (20) and the dashed line (9) is exactly the factor of two due to the extra degree of free- 
dom derived in Section 3.1. At small radii (r < 5 kpc), the solutions diverge: the numerical 
solution contains rings of the same radius but different tilts while the simulation bins the 
particles by radius and therefore averages particles in the same shell together; more impor- 
tantly, the assumption that dr/dt is small is violated. However, in the innermost region of a 
halo the symmetric bulge will dominate the potential, and therefore the torque power law is 
likely invalid in any case. For these reasons, we ignore the innermost 2 kpc in our analysis 
and consider equation (20) to be a good description of the disk everywhere else. 



4. Self-gravitating Disks and Warps 

Figure 4 suggests that disks should have continuously curved warps in their inner regions, 
eventually becoming flat in the outer regions where t t n t is small. Real warped galaxies 
appear quite the opposite; they are solid out to approximately the Holmberg radius (half 
of the diameter of the I = 26.5 /i pg isophote) and warped beyond that (Briggs 1990). The 
most important difference between real galactic disks and those used in section 3 is that the 
self-gravity of real galactic disks is important; galaxies have mass. 

Although the dark matter halo must dominate the potentials of disk galaxies at large 
radii, there is considerable evidence that the self-gravity of the disk is important to its 
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dynamics over much of its area. Broeils & Courteau (1997) show that "maximal disk" 
models, where the mass-to-light ratios of the disks (M/L) are as large as is consistent with 
the rotation curves, provide good fits to most rotation curves. Sackett (1997) demonstrates 
that the Milky Way is also consistent with having a maximal disk. Although Courteau & Rix 
(1999) have used the residuals of the Tully-Fisher relation (Tully & Fisher 1977) to estimate 
that disks provide only 40% of the dynamical mass at 2.2 exponential scale lengths, and de 
Blok & McGaugh (1997) have demonstrated that the rotation curves of low surface brightness 
galaxies imply that they are substantially sub-maximal, the prevalent phenomena of spiral 
arms and bars, which are disk instabilities, require that in many disks, the self-gravity of the 
disk must be important to its dynamics (Athanassoula et al 1987). 

The self-gravity of the disk acts to keep the disk flat. Ostriker & Binney (1989) examined 
the effect of a slewing disk potential on a set of self-gravitating rings and found that regions 
of high surface density react like a solid body. Since the surface density of the disk is highest 
in the central regions, the central parts of the disk will be kept locally flat, resulting in 
disks that closer resemble observed warped galaxies. Using numerical N-body simulations, 
we investigate this effect. 

Lopez-Corredoira et al (2002a) take this into account by developing a form for the 
internal torque Tint- Here we do not use their r- mt , but rather integrate the orbits numerically. 
While this does not provide us with an explicit set of differential equations for 8(r,t), it 
frees us from such assumptions as the constancy of u(r) (which we saw was important in 
Section 3.1), the existence of an equilibrium configuration, and the lack of a dark matter 
halo. 



4.1. The simulations 

To evaluate how the mass of the disk affects the results of section 3, we performed 
numerical N-body simulations of massive Milky Way type galactic disks. Equilibrium disk 
models of disk mass 1 x 1O 1O M , 3 x 10 10 M Q , and 5.6 x 1O 1O M , scale length r d = 3.5 kpc 
and vertical scale height h z = 325 pc in a static spherically-symmetric NFW halo potential 
(Navarro et al 1997) with concentration parameter C200 = 15 and virial velocity V200 = 
175 km s -1 were constructed using the method of Hernquist (1993) and then allowed to relax 
under the force of gravity until they appeared to be in equilibrium. The halo masses were 
varied such that the mass of the disk plus halo was 4.1 x 10 11 M in each case. 

The models were evolved using the GRAPESPH code (Steinmetz 1996), where inter- 
particle gravitational forces are computed using direct summation with GRAPE-3 hardware 
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(Okumura et al 1993). A plummer softening of 0.3 kpc has been used. The models were 
evolved for 2 Gyr, which took 5000-7000 timesteps depending on the model. The torque 
force was applied as an external force as described in section 3.2 with r = 10~ 30 s~ 2 at 
r = 1 kpc. Simulations were performed with torque slopes of (3 — —2.0, —2.5, and —3.0. 

The disks contained 16384 particles. The 3 x 10 10 M simulations were also performed 
with 32767 particles to see if the resolution was sufficient. The results for the high resolution 
simulations were identical to those for the lower resolution simulations to well within the 
angular errors computed in section 4.3 except for in the spherical shells which contained very 
few particles and had been neglected from the analysis on those grounds. 

4.2. Warped disks 

Figure 5 shows the simulation of a 3 x 10 10 M disk with 32767 particles subject to a 
j3 = —2.5 torque for 1 Gyr. The main plane of the disk is flat and clearly tilted toward the 
positive x-axis, as expected from Section 2.1.2. Beyond 10 kpc, the disk no longer remains 
flat but warps back toward the original plane. The particles which appear to be filling in 
the area between the main disk and the warp are projections and are actually in front of or 
behind the galaxy at large radii. 

We binned the simulation particles into spherical shells of thickness 2 kpc and calculated 
their tilts as in section 3.2. Figure 6 shows a plot of the radial profile of the disk shown 
in Figure 5, along with the massless simulation and analytic prediction (20) from Figure 4, 
which correspond to a massless disk in the same torque. There are three important regimes: 
the central flat disk, the warp, and the effectively massless outer region. The central 11 kpc 
is all tilted 17° from the original plane (we ignore the innermost point because the torque 
is capped inside tq), and corresponds to the visually flat part of the disk. The effect of the 
disk's self-gravity is to prevent the inner regions of the disk from tilting as much as they 
would otherwise, while pulling the outer regions of the disk to larger tilts. At the warp radius 
r w =ll kpc, the disk is no longer flat, and the disk warps back toward the original plane. 
Finally, from 17 kpc to 23 kpc (the extent of the disk), it follows the analytic prediction, 
acting like a massless disk. This result agrees qualitatively with that of Ostriker & Binney 
(1989), who found that regions of high surface density remain flat when torqued, and the 
position of a warp is determined in part by a drop in surface density. 



-16- 




Fig. 5. — x-z projection of a simulated disk galaxy with mass 3 x 10 M© after 1 Gyr in a 
torque of r = 10 -30 s~ 2 at r = 1 kpc and (3 = —2.5. 
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Fig. 6. — Radial dependence of the tilt of the disk in Figure 5, computed in spherical shells 
2 kpc wide. The solid line is the analytic prediction (20) and the solid diamonds are the 
massless N-body simulation, as in Figure 4. 
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4.3. Warp evolution 

We now investigate how the disk evolves over time under the influence of the torque. 
The general form of the warp is as above, but as the disk evolves under the influence of the 
torque, the warp radius moves out through the disk at a rate that depends on the mass of 
the disk. At early times, typically before 200 Myr, the solid region is not well developed and 
lies below the massless prediction. 

We developed an algorithm to automatically detect the warp radius. The particle posi- 
tions were stored every 40 timesteps. For each of these outputs, the tilt angles were calculated 
in spherical shells as in Figure 6. We did a bootstrap analysis to estimate the error in these 
angles: for each spherical shell, we drew 10 random sets from the particles in that shell, cal- 
culated the tilt of each bootstrap set, and then used the standard deviation of those angles 
as an estimate of the error in the tilt of that bin. The error depended on the number of 
particles in the shell, but was typically less than 1°. 

Shells with errors greater than 0.5° were not analyzed, as these generally had few par- 
ticles and were dominated by numerical noise. For each remaining shell with radius r and 
tilt # S i m , we calculated the deviation A6 between the simulation tilt # S i m and the analytic 
prediction 9 pre< i(r) from equation (20). The average radius r of the radial bin with the 
highest A8 was then considered to be the warp radius r w . This is the radius at which the 
difference between the tilt of the massive disk and the tilt of an equivalent massless disk 
is maximized. For example, in Figure 6, the highest A6 occurs at r = 11 kpc where the 
data point # s j m = 17° and the predicted line #p rec [ = 7°. This corresponds exactly to the 
radial bin at the end of the plateau in Figure 6. Examining a few randomly chosen outputs 
from each simulation showed that in each case this automated r w agreed with our intuition, 
except at early times when there was no solid disk and the automated r w was dominated by 
noise. 

Figures 7, 8, and 9 show the evolution of this warp radius for the set of simulations 
with differing torque slopes (5. The squares, circles, and triangles correspond to simulations 
with disk masses of 1, 3, and 5.6 x 10 10 M Q respectively. The warp radii appear "quantized" 
because of the 2 kpc- wide radial bins used to calculate the tilt (see, e.g. Figure 6). One 
point was taken every 10 s years, linearly interpolating r w from the neighbouring simulation 
outputs. The early times before the warp has developed show up as noise in the upper-left 
region of some of the figures, particularly Figure 9. To quantify the growth of the warp, 
straight lines were fit through the profiles. While one might think that the lines ought to 
be constrained to pass through the origin, since there is no warp at t — 0, the data do not 
support this because the warp instability requires that the surface density of the disk be 
below a certain critical surface density, as we demonstrate later. We examined Figures 7-9 
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Fig. 7. — Radial bins at which the disks exhibit their warps every 10 8 yr as the simulations 
evolve in time. Because the particles were binned into spherical shells 2 kpc wide, the 
warp radii appear quantized. The lines are linear least square fits. This figure shows the 
simulations with f3 = —2.0 torques and disk masses 1.0 x 10 10 M Q , 3.0 x 10 10 M , and 
5.6 x 10 10 M & represented by the squares, circles, and triangles respectively. 
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Fig. 8. — As in Figure 7 but for simulations with (3 = —2.5 torques. 
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Fig. 9. — As in Figure 7 but for simulations with (5 = —3.0 torques. The warp takes time to 
develop, resulting in the noise apparent early in the simulations, which were not included in 
the straight line fits. 
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by eye and excluded the obviously noisy early times from the fits; in practice, this meant 
excluding t < 5 x 10 s yr for the (3 = -3.0, M d = 1.0 x 10 10 M simulation and t < 2 x 10 8 yr 
for the f3 = —3.0, = 3.0 x 10 10 M simulation. We ended the simulations after 2 Gyr, the 
time when solid region reached the edge of the particle disk in some simulations, at which 
point the algorithm to find r w just detects the edge of the disk. The fits are shown as the 
solid, dashed, and dotted lines for the simulations with disk masses of 1, 3, and 5.6 x 1O 1O M 
respectively. 

The slope of the linear fit measures how quickly the warp grows. This could depend 
both on the mass of the disk, which affects the local surface density at the warp radius, and 
the slope of the torque, which determines the strength of the torque at the warp radius. 
Figure 10 shows the rate that r w grows for each simulation. There is an apparent trend that 
the warps in higher mass disks move through the disk faster than in lower mass disks, but 
the scaling does not appear simple. There is some indication from Figure 10 that the warps 
from shallower torque laws (which have stronger torques at the warp radii) grow faster than 
those with steep torques, but the situation is not clear. 

The effect of the disk mass on the rate of the warp growth implies that the self-gravity 
of the disk is important to the formation of the warp. In a massive disk, particles at different 
radii are coupled to each other gravitationally and act to keep the disk locally flat. This 
suggests that the crucial factor that determines where the warp develops is the local surface 
density of the disk. Ostriker & Binney (1989) examined the effect of a slewing disk potential 
on a set of self-gravitating rings and found qualitatively that regions of high surface density 
react like a solid body, but that warps can occur where the surface density is lower. Hofner 
& Sparke (1994) noted that in most cases the group speed of bending waves in a disk of 
surface density E(r) and angular rotation velocity uj{r) is 

7rGS(r) 

and therefore the time for a warp to settle at a given radius would be inversely proportional 
to the surface density at that radius. 

To test this, we translate the warp radii of the simulations into local surface densities 
using 

for a disk of mass M d and exponential scale length r e . Figure 11 shows the surface density 
at the warp radius as a function of time for the three disks of different mass in a j3 = —2.5 
torque (i.e. for the same simulations shown in Figure 8), while Figure 12 shows it for a 
3.0 x 10 10 M disk in torques with varying (3. The surface density at the warp radius falls 
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Fig. 10. — Rate that the warp moves out through the disk (i.e. the slopes of the linear fits 
in Figures 7-9) as a function of disk mass, with different torque slopes (3 denoted by the 
different symbols. The effect of the slope on the rate of warp growth is ambiguous, but there 
is a clear trend that more massive disks have faster-growing warps. 
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Fig. 11. — Surface densities at the warp radii as the simulations evolve. As the warp moves 
out through the disk, the local surface density at that radius falls according to equation (23). 
The squares, circles, and triangles show simulations with disk masses of 1, 3, and 5.6 x 1O 1O M 
respectively in a (3 = —2.5 torque. The line is the fit (24). Compare this with Figure 8, 
which shows the radial evolution of the warp for the same set of simulations. The local 
surface density at the warp is similar for all models at all times. 
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Fig. 12. — As in Figure 11, but plotting 3.0 x 1O 1O M disk simulations with different torque 
slopes (3. 
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as the warp moves out through the disk. The local surface densities at the warp at a given 
time are quite similar for all disk masses, and have even smaller scatter for different torque 
laws. It appears that the local surface density is the important parameter for determining 
the warp radius at a given time. 

The evolution of the warp is well described by a decaying exponential 

E(r^) = £ Wo e-^ (24) 

with T, Wo = 70 M Q pc~ 2 and t w = 480 Myr, shown as the straight line in Figures 11 and 12. 
For such a result to hold, the timescale t w should depend only on global properties of the 
galaxy independent of disk mass or torque; in particular, it can only depend on E wo , the virial 
velocity of the galaxy t> 2 oo = 175 kms^" 1 , the exponential scale length of the disk = 3.5kpc, 
and/or the vertical scale height h z = 325 pc. An interesting timescale that matches this is 
the characteristic timescale for bending waves (Hofner & Sparke 1994) one disk scale length 
away from E wo : 

t g = -= V200 (25) 
c„ 7iGYj W0 e 



which is 490 Myr in these cases. 



H WQ is the extrapolation of equation (24) to t — 0, and is the critical surface density 
above which the disk does not develop a warp. At higher surface densities, the self-gravity 
of the disk is always sufficient to keep the disk flat. It is interesting that the extrapolation to 
t — 0, when there is no torque, gives a finite well-defined value for the warp surface density. 
This is why the linear r w (t) fits do not go through the origin and suggests that at this 
surface density the disk is marginally unstable to warping. We do not expect to see warps 
occurring at surface densities higher than 70 M pc~ 2 . In the Milky Way, the warp begins 
at or slightly beyond the solar circle (Binney 1992), while recent Hipparcos determinations 
of the Milky Way surface density in the solar neighbourhood give E = 40 M pc -2 (Creze 
et al 1998). It is encouraging that this is less than E wo . In external galaxies, warps begin 
between 25 and 26.5 mag arcsec -2 (Briggs 1990), which for a mass to light ratio in B of 
1M & /L & gives surface densities of between 1.8-7.0 M Q pc -2 , ignoring projection effects and 
extinction. This is also consistent with a picture in which warps only can occur at surface 
densities below E„ 



J WQ- 



5. Summary 



We studied how a typical galactic disk reacts to torques expected from lying in mis- 
aligned dark matter halos. Our main findings are the following: 
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• Cosmological N-body simulations suggest that galactic disks will be misaligned with 
the mass distribution of the dark matter halo in which they are embedded. Because 
of this misalignment, the halo exerts a perpendicular gravitational torque on the disk 
of the form 



^ I 3 



1 kpc, 

with typical values of t = lCT 30 s~ 2 and j3 = —2.5. This gravitational torque is strong 
enough to have a significant effect on the entire disk within a Hubble time. 

The timescale for an orbit to tilt in response to this torque rises with radius, so inner 
portions of the disk will realign themselves first, resulting in a warped disk. For a 
massless disk of stars in circular orbits, a good estimate of how tilted the disk is at a 
radius r after a time t is 

1 l + e */*°W 



6(r,t) = 2cos _1 
where the timescale t (r) is defined as 

to(r) 



2v r (r) ( r x 



Tor \r 0/ 

and is twice as long as in the case where the orbital radii are not allowed to vary. 

• Massive disks depart from this due to the self-gravity of the inner portions of the disk. 
The disk is kept flatter where the local surface density is high than for a massless disk. 
The radius inside which the disk is flat (which corresponds to the radius where the 
warp would be considered to start if the disk were observed) grows with time. The 
warp grows at a rate of between 3 and 10 kpc Gyr -1 . The rate depends on the mass 
of the disk, but is relatively insensitive to the torque parameters. 

• More massive disks have faster-growing warps, since only at large radius is the surface 
density sufficiently low that its self-gravity cannot maintain the flat disk. The surface 
density at which the warp occurs is well-described by 

where S wo = 70 M Q pc~ 2 and t w = 480 Myr, independent of disk mass and torque 
slope. is the maximum surface density at which a warp can occur, and marks the 
point where a disk is marginally unstable to warping. t w appears to coincide with the 
timescale of bending waves one disk scale length away from Y1 WQ : 
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In the future, we plan to analyze high resolution cosmological simulations to determine the 
torques due to realistic dark matter halos. This will enable us to distinguish the torque due 
to the misaligned disk, that due to substructure in the halo, and that due to the external 
tidal field by modeling individual halos. A further direction that we have not yet explored 
is sampling cosmological simulations at regular intervals up to the present day to replace 
the static torque force used in this work with a torquing field that changes over time in a 
cosmologically-motivated manner. 

Several figures were produced using the WIP package (Morgan 1995). This work was 
supported by NSF grants 9807151 and PHY99-0749, and NSERC grant PGSB-233028-2000. 
MS is an Alfred P. Sloan Fellow and a David and Lucile Packard Fellow. 
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